Solution of the tunneling-percolation problem in the nanocomposite regime 
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We noted that the tunnehng-percolation framework is quite well understood at the extreme cases 
of percolation-like and hopping-like behaviors but that the intermediate regime has not been previ- 
ously discussed, in spite of its relevance to the intensively studied electrical properties of nanocom- 
posites. Following that we study here the conductivity of dispersions of particle fillers inside an insu- 
lating matrix by taking into account explicitly the filler particle shapes and the inter-particle electron 
tunneling process. We show that the main features of the filler dependencies of the nanocomposite 
conductivity can be reproduced without introducing any a priori imposed cut-off in the inter-particle 
conductances, as usually done in the percolation-like interpretation of these systems. Furthermore, 
we demonstrate that our numerical results are fully reproduced by the critical path method, which is 
generalized here in order to include the particle filler shapes. By exploiting this method, we provide 
simple analytical formulas for the composite conductivity valid for many regimes of interest. The 
validity of our formulation is assessed by reinterpreting existing experimental results on nanotube, 
nanofiber, nanosheet and nanosphere composites and by extracting the characteristic tunneling de- 
cay length, which is found to be within the expected range of its values. These results are concluded 
then to be not only useful for the understanding of the intermediate regime but also for tailoring 
the electrical properties of nanocomposites. 
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I. INTRODUCTION 

The inclusion of nanomctric conductive fillers such as 
carbon nanotubesi, nanofibcrsS., and grapheno^ into in- 
sulating matrices allows to obtain electrically conductive 
nanocomposites with unique properties which are widely 
investigated and have several technological applications 
ranging from antistatic coatings to printable electronics^. 
A central challenge in this domain is to create compos- 
ites with an overall conductivity a that can be controlled 
by the volume fraction 0, the shape of the conducting 
fillers, their dispersion in the insulating matrix, and the 
local inter-particle electrical connectedness. Understand- 
ing how these local properties affect the composite con- 
ductivity is therefore the ultimate goal of any theoretical 
investigation of such composites. 

A common feature of most random insulator-conductor 
mixtures is the sharp increase of a once a critical volume 
fraction (j)c of the conductive phase is reached. This tran- 
sition is generally interpreted in the framework of perco- 
lation theory^"— and associated with the formation of a 
cluster of electrically connected filler particles that spans 
the entire sample. The further increase of a for (f) > (f>c 
is likewise understood as the growing of such a cluster. 
In the vicinity of (f)c, this picture implies a power-law 
behavior of the conductivity of the form 



(1) 



where f is a critical exponent. Values of t extracted from 
experiments range from its expected universal value for 
three-dimensional percolating systems, i ~ 2, up to i ~ 
10, with little or no correlation to the critical volume 
fraction (per or the shape of the conducting fillers^. 



In the dielectric regime of a system of nanometric 
conducting particles embedded in a continuous insulat- 
ing matrix, as is the case for conductor-polymer nano- 
compositesr^"— the particles do not physically touch 
each other, and the electrical connectedness is estab- 
lished through tunneling between the conducting filler 
particles. In this situation, the basic assumptions of 
percolation theory are, a priori, at odds with the inter- 
particle tunneling mechanismiiiii Indeed, while percola- 
tion requires the introduction of some sharp cut-off in the 
inter-particle conductances, i.e., the particles are either 
connected (with given non-zero inter-particle conduc- 
tances) or disconnccted^i^ the tunneling between par- 
ticles is a continuous function of inter-particle distances. 
Hence, the resulting tunneling conductance, which de- 
cays exponentially with these distances, does not imply 
any sharp cut-off or threshold. 

Quite surprisingly, this fundamental incompatibility 
has hardly been discussed in the literature?^ and ba- 
sically all the measured conductivity dependencies on 
the fractional volume content of the conducting phase, 
ct((/)), have been interpreted in terms of Eq. ([T]) assum- 
ing the "classical" percolation behavior In this article, 
we show instead that the inter-particle tunneling explains 
well all the main features of cr(0) of nanocomposites with- 
out imposing any a priori cut-off, and that it provides a 
much superior description of cr{(f)) than the "classical" 
percolation formula ([T|). 

In order to specify our line of reasoning and to bet- 
ter appreciate the above mentioned incompatibility, it is 
instructive to consider a system of particle dispersed in 
an insulating continuum with a tunneling conductance 
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between two of them, i and j, given by: 

26, 



goexp 



(2) 



where go is a constant, ^ is the characteristic tunneUng 
length, and 6ij is the minimal distance between the two 
particle surfaces. For spheres of diameter D, 5ij ~ rij~D 
where is the center-to-center distance. There are two 
extreme cases for which the resulting composite conduc- 
tivity has qualitatively different behaviors which can be 
easily described. In the first case the particles are so 
large that ^/D 0. It becomes then clear from Eq.® 
that the conductance between two particles is non-zero 
only when they essentially touch each other. Hence, re- 
moving particles from the random closed packed limit is 
equivalent to remove tunneling bonds from the system, 
in analogy to sites removal in a site percolation problem 
in the lattice.— The system conductivity will have then 
a percolation-like behavior as in Eq. ^ with t ~ 2 and 
0c being the corresponding percolation thresholdJ^ The 
other extreme case is that of sites (D/^ — ^ 0) randomly 
dispersed in the continuum. In this situation, a variation 
of the site density p does not change the connectivity 
between the particles and its only role is to vary the dis- 
tances 5ij = Tij between the sites i^^'^^ The correspond- 
ing fj behavior was solved by using the critical path (CP) 
methodic in the context of hopping in amorphous semi- 
conductors yielding a oc cxp[— 1.75/(^p^/'^)]ii^iiS For suf- 
ficiently dilute system of impenetrable spheres this rela- 
tion can be generalized to cr oc exp[— 1.41D/ (^0^/^)] It 
is obvious then from the above discussion that the second 
case is the low density limit of the first one, but it turns 
out that the variation of (t(0) between the two types of 
situations, which is definitely pertinent to nanocompos- 
ites, has not been studied thus far. 

Following the above considerations we turned to study 
here the cr((/)) dependencies by extending the low-density 
(hopping-like) approach to higher densities than those 
used previously. ^^'^^1-'^^ Specifically, we shall present nu- 
merical results obtained by using the global tunneling 
network (GTN) model, where the conducting fillers form 
a network of globally connected sites via tunneling pro- 
cesses. This model has already been introduced in 
Ref.[20j for the case of impenetrable spheres, but here 
we shall generalize it in order to describe also anisotropic 
fillers such as rod-like and plate-like particles, as to ap- 
ply to cases of recent interest (i.e., nanotube, nanofiber, 
nanosheet, and graphene composites). In particular, the 
large amount of published experimental data on these 
systems allows us to test the theory and to extract the 
values of microscopic parameters directly from macro- 
scopic data on the electrical conductivity. 

The structure of the paper is as follows. In Sec. [Tl] 
we describe how we generate particle dispersions and in 
Sec. mil we calculate numerically the composite conduc- 
tivities within the GTN model and compare them with 
the conductivities obtained by the CP approximation. In 
Sec. IIVI we present our results on the critical tunneling 



distance which are used in Sec.|V]to obtain analytical for- 
mulas for the composite conductivity. These are applied 
in Sec. I VII to several published data on nanocomposites 
to extract the tunneling distance. Section IVlIl is devoted 
to discussions and conclusions. 



II. SAMPLE GENERATION 

In modeling the conductor-insulator composite mor- 
phology, we treat the conducting fillers as identical im- 
penetrable objects dispersed in a continuous insulating 
medium, with no interactions between the conducting 
and insulating phases. As pointed out above, in order 
to relate to systems of recent interest we describe filler 
particle shapes that vary from rod-like (nanotubes) to 
plate-like (graphene). This is done by employing impen- 
etrable spheroids (ellipsoids of revolution) ranging from 
the extreme prolate {a/b ^ 1) to the extreme oblate 
limit (a/6 ^ 1), where a and h are the spheroid polar 
and equatorial semi-axes respectively. 

We generate dispersions of non-overlapping spheroids 
by using an extended version of a previously described 
algorithm^! which allows to add spheroids into a cubic 
cell with periodic boundary conditions through random 
sequential addition (RSA) 3^ Since the configurations ob- 
tained via RSA are non-equilibrium ones,—— the RSA 
dispersions were relaxed via Monte Carlo (MC) runs, 
where for each spheroid a random displacement of its cen- 
ter and a random rotation of its axes2^ were attempted, 
being accepted only if they did not give rise to an over- 
lap with any of its neighbors. Equilibrium was consid- 
ered attained when the ratio between the number of ac- 
cepted trial moves versus the number of rejected ones 
had stabilized. Furthermore, to obtain densities beyond 
the ones obtainable with RSA, a high density generation 
procedure22i2£ was implemented where in combination 
with MC displacements the particles were also inflated. 
The isotropy of the distributions was monitored b y u sing 
the nematic order parameter as described in Ref.|27|. 

Figure [1] shows examples of the so-generated distribu- 
tions for spheroids with different aspect-ratios a/b and 
volume fractions (p = V where V = 47ra6^/3 is the vol- 
ume of a single spheroid and p is the particle number 
density. 



III. DETERMINATION OF THE COMPOSITE 
CONDUCTIVITY BY THE GTN AND CP 
METHODS 

In considering the overall conductivity arising in such 
composites, we attributed to each spheroid pair the tun- 
neling conductance given in Eq. ^ where, now, for 
a/b ^ 1 the inter-particle distance 5ij depends also on 
the relative orientation of the spheroids. The 5ij val- 
ues were obtained here from the numerical procedure de- 
scribed in Ref. 2l|. On the other hand, in writing Eq. ^ 



3 




FIG. 1: (Color online) Examples of distributions of impenetrable spheres and spheroids of different aspect-ratios a/b and volume 
fraction generated by the algorithms used in the present work. 



we neglect any energy difference between spheroidal par- 
ticles and disregard activation energies since, in gen- 
eral, these contributions can be ignored at relatively high 
temperatures,— which is the case of interest here. For 
the specific case of extreme prolate objects [a/h 3> 1) the 
regime of validity of this approximation has been studied 
in Ref.@. 

The full set of bond conductances given by Eq. ^ was 
mapped as a resistor network with go = 1 and the overall 
conductivity was calculated through numerical decima- 
tion of the resistor network.i^^*^ To reduce computational 
times of the decimation procedure to manageable limits, 
an artificial maximum distance was introduced in order 
to reject negligibly small bond conductances. It is impor- 
tant to note that this artifice is not in conflict with the 
rationale of the GTN model, since the cutoff it implies 



neglects conductances which are completely irrelevant for 
the global system conductivity. We chose the maximum 
distance to be generally fixed and equal to four times 
the spheroid major axis (i.e. a in the prolate case and 
h in the oblate case), which is equivalent to reject inter- 
particle conductances below e"^" for ^/D = 1/15 case 
(and considerably less for smaller ^ values). However, 
for the high aspect-ratios and high densities the distance 
had to be reduced. Moreover, since the maximum dis- 
tance implies in turn an artificial geometrical percolation 
threshold of the system, for the high aspect-ratios, at 
low volume fractions the distance had to be increased to 
avoid this effect. By comparing the results with the ones 
obtained with significantly larger maximum distances we 
verified that the effect is undetectable. 

In Fig. HJa) we show the so-obtained conductivity a 
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FIG. 2: (Color online) The results of our GTN and CP calculations, (a): Volume fraction 4> dependence of the tunneling 
conductivity a for a system of aspect-ratio a/h = 10 hard prolate spheroids with different characteristic tunneling distances 
^/Z) with D = 2a. Results from Eq. ^ (with cto = 0.179) are displayed by dotted lines, (b): Tunneling conductivity in a 
system of hard spheroids with different aspect-ratios a/h and (^/D = 0.01, with D — 2max(a,&). Dotted lines: results from 
Eq. (© with (70 = 0.124 for a/b = 2, ao = 0.099 for a/b = 1/2, ao = 0.351 for a/b = 1/10, and ctq = 0.115 for a/b = 1. 



values (symbols) as a function of the volume fraction 
(j) of prolate spheroids with aspect-ratio a/b = 10 and 
different values of ^/D, where D = 2max(a, 6). Each 
symbol is the outcome of = 200 realizations of a 
system of Np ^ 1000 spheroids. The logarithm average 
of the results was considered since, due to the exponen- 
tial dependence of Eq. the distribution of the com- 
puted conductivities was approximately of the log- normal 
form^i The strong reduction of a for decreasing (j) shown 
in Fig. [IJa) is a direct consequence of the fact that as 4> is 
reduced, the inter-particle distances get larger, leading in 
turn to a reduction of the local tunneling conductances 
[Eq. ([2])]. In fact, as shown in Fig. [U^b), this reduction 
depends strongly on the shape of the conducting fillers. 
Specifically, as the shape anisotropy of the particles is en- 
hanced, the composite conductivity drops for much lower 
values of (j) for a fixed ^. 

Having the above result we turn now to show that 
the strong dependence of cr((/)) on a/b and in Fig. [2] 
can be reproduced by CP methodi^"— when applied to 
our system of impenetrable spheroids. For the tunneling 
conductances of Eq. this method amounts to keep 
only the subset of conductances gtj having dij < Sc, 
where Sc, which defines the characteristic conductance 
5c — 5o exp(— 25c/C); is the largest among the Sij dis- 
tances, such that the so-defined subnetwork forms a con- 
ducting cluster that span the sample. Next, by assigning 
gc to all the (larger) conductances of the subnetwork, a 
CP approximation for a is 

25c{(f>,a,b)] 
cr ~ Cto exp , (3) 

where ctq is a pre- factor proportional to go- The signif- 
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FIG. 3: (Color online) Schematic illustration of the tunneling 
conductivity crossover for the cases a/b — 1, a/b = 2, and 
a/b = 10. 



icancc of Eq. ([3]) is that it reduces the conductivity of 
a distribution of hard objects that are electrically con- 
nected by tunneling to the computation of the geometri- 
cal "critical" distance Sc- In practice, Sc can be obtained 
by coating each impenetrable spheroid with a penetrable 
shell of constant thickness 6/2, and by considering two 
spheroids as connected if their shells overlap. Sc is then 
the minimum value of i5 such that, for a given 0, a cluster 
of connected spheroids spans the sample. 

To extract Sc we follow the route outlined in Ref. [2l| 
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with the extended distribution generation algorithm de- 
scribed in Sec. [Ill Specificahy, we calculated the span- 
ning probability as a function of cj) for fixed a/b and 6c 
by recording the frequency of appearance of a percolat- 
ing cluster over a given number of realizations -/V^. The 
realization number varied from Nfi = 40 for the smallest 
values of 6c up to Nfj = 500 for the largest ones. Each re- 
alization involved distributions of Np ^ 2000 spheroids, 
while for high aspect-ratio prolate spheroids this number 
increased to Np ~ 8000 in order to be able to main- 
tain the periodic boundary conditions on the simulation 
cell. Relative errors on 6c were in the range of a few per 
thousand. 

Results of the CP approximation are reported in Fig. [5] 
by dotted lines. The agreement with the full numerical 
decimation of the resistor network is excellent for all val- 
ues oi a/b and £,/D considered. This observation is quite 
important since it shows that the CP method is valid also 
beyond the low-density regime, for which the conducting 
fillers are effectively point particles, and that it can be 
successfully used for systems of particles with impene- 
trable volumes. Besides the clear practical advantage of 
evaluating a via the geometrical quantity 6c instead of 
solving the whole resistor network, the CP approxima- 
tion is found then, as we shall see in the next section, to 
allow the full understanding of the filler dependencies of 
a and to identify asymptotic formulas for many regimes 
of interest. 

Before turning to the analysis of the next section, it 
is important at this point to discuss the following is- 
sue. As shown in Fig. [21 the GTN scenario predicts, 
in principle, an indefinite drop of cr as ^ — >■ be- 
cause, by construction, there is not an imposed cut- 
off in the inter-particle conductances. However, in real 
composites, either the lowest measurable conductivity is 
limited by the experimental set-up,— or it is given by 
the intrinsic conductivity CTm of the insulating matrix, 
which prevents an indefinite drop of a. For example, in 
polymer-based composites CTm falls typically in the range 
of (Tm — 10^^'^ ^ lO^^'* S/cm, and it originates from ionic 
impurities or displacement currents.— Since the contri- 
butions from the polymer and the inter-particle tunneling 
come from independent current paths, the total conduc- 
tivity (given by the polymer and the inter-particle tun- 
neling) is then simply crtot = fm + cr.— As illustrated 
in Fig. [31 where trtot is plotted for a/b = 1, 2, and 10 
and for am/co = 10~^^, the 0-dependence of crtot is char- 
acterized by a cross-over concentration (j)c below which 
CTtot — fm- As seen in this figure, fillers with larger 
shape- anisotropy entail lower values of (/)c, consistently 
with what is commonly observedjii^"— We have there- 
fore that the main features of nanocompositcs (drop of a 
for decreasing (f>, enhancement of a at fixed (j) for larger 
particle anisotropy, and a characteristic (j)c below which 
the conductivity matches that of the insulating phase) 
can be obtained without invoking any microscopic cut- 
off, leading therefore to a radical change of perspective 
from the classical percolation picture. In particular, in 



the present context, the conductor- insulator transition is 
no longer described as a true percolation transition (char- 
acterized by a critical behavior of a in the vicinity of a 
definite percolation threshold, i.e., Eq.(|T])), but rather as 
a cross-over between the inter-particle tunneling conduc- 
tivity and the insulating matrix conductivity. 

IV. CP DETERMINATION OF THE CRITICAL 

DISTANCE 5c FOR SPHEROIDS 

The importance of the CP approximation for the un- 
derstanding of the filler dependencies of a is underscored 
by the fact that, as discussed below, for sufficiently elon- 
gated prolate and for sufficiently flat oblate spheroids, as 
well as for spheres, simple relations exist that allow to 
estimate the value of 6c with good accuracy. In virtue of 
Eq. ([3]) this means that we can formulate explicit rela- 
tions between a and the shapes and concentration of the 
conducting fillers. 

A. prolate spheroids 

Let us start with prolate {a/b > 1) spheroids. In 
Fig. [H^a) we present the calculated values of 6c/ D as a 
function of the volume fraction for spheres (a/b = 1, 
together with the results of Ref. [s^l) and for a/b = 2, 
10, 20, and 100. In the log-log plot of Fig. \M})) the 
same data are displayed with 6c/ D multiplied by the ra- 
tio V^phcrc/V^ = (a/&)^, where T^phcrc = ttD^/6 is the 
volume of a sphere with diameter equal to the major 
axis of the prolate spheroid and V = Airab'^/S is the 
volume of the spheroid itself. For comparison, we also 
plot in Fig. HKb) the results for impenetrable sphero- 
cylinders of Refs. [27ll38| . These are formed by cylin- 
ders of radius R and length L, capped by hemispheres 
of radius i?, so that a = R + L/2 and b = R, and 

V, phcrc/V = {a/bf/\{2>/2){a/b) - 2] ~ (2/3)(a/6)2 for 
a/b ;3> 1. As it is apparent, for sufficiently large values of 
a/b the simple re-scaling transformation collapses both 
spheroids and spherocylinders data into a single curve. 
This holds true as long as the aspect-ratio of the spheroid 
plus the penetrable shell (a 4- (5c/2)/(6 -f- (5c/2) is larger 
than about 5. In addition, for (\) < 0.03 the collapsed data 
are well approximated by ^cV^phoro/^^/i^ = 0.4/0 [dashed 
line in Fig.[4l^b)], leading to the following asymptotic for- 
mula: 

W^.^. (4) 

where 7 = 0.4 for spheroids and 7 = 0.6 for spherocylin- 
ders. Equation (|4|) is fully consistent with the scaling 
law of Ref. [S^ that was obtained from the second- virial 
approximation for semi-penetrable spherocylinders, and 
it can be understood from simple excluded volume ef- 
fects. Indeed, in the asymptotic regime a/b ^ 1 and for 
6c/ a ^ 1, the filler density p (such that a percolating 
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FIG. 4: (Color online) (a): the 5c/ D dependence on the volume fraction for im pen etrable prolate spheroids with a/h = 1, 2, 
10, 20, and 100. For a/h — 1 our results are plotted together with those of Ref. [33 . The solid line is E g. jO jl. (b): re-scaled 
critical distances versus for prolate spheroids as well as for the impenetrable spherocylinders of Refs. [27ll3q . The dashed 
line follows Eq. Q and the solid line follows Eq. (O 



cluster of connected semi-penetrable spheroids with pen- 
etrable shell 5c is formed) is given by p = l/AVcxc-— '^'^ 
Here, AVcxc is the excluded volume of a randomly ori- 
ented scmi-penctrable object minus the excluded volume 
of the impenetrable object. As shown in the Appendix, 
for both spheroids and spherocylindcr particles this be- 
comes AVoxc — 2710? 5c, leading therefore to Eq. ^ with 
7 = 1/3 for spheroids and 7 = 1/2 for spherocylinders4^ 
It is interesting to notice that in Fig. |4l[b) the re-scaled 
data for > 0.03 deviate from Eq. (jj]) but still follow a 
common curve. We have found that this common trend 
is well fitted by an empirical generalization of Eq. 



5c/D 



lib/a) 
Ml + 8' 



(5) 



which applies to all values of (/> provided that (a 
5c/2)/{b + 5c/2) > 5 [solid lines in Fig.Hb)]. 



B. oblate spheroids 

Let us now turn to the case of oblate spheroids. The 
numerical results for 5c as a function of the volume frac- 
tion (j) are displayed in Fig. [Sja) for a/b = 1, 1/2, 1/10, 
1/100, and 1/200. Now, as opposed to prolate fillers, 
almost all of the experimental results on nanocompos- 
ites, such as graphenc;^ that contain oblate filler with 
high shape anisotropy are at volume fractions for which 
a corresponding hard spheroid fiuid at equilibrium would 
already be in the nematic phase. For oblate spheroids 
with a/b = 1/10 the isotropic- nematic (I-N) transition is 
at 0i_N ~ 0.185)^ while for lower a/b values the tran- 
sition may be estimated from the results on infinitely 



thin hard disks44 0i_n ^ 0.0193 for a/b = 1/100 and 
^i_N ^ 0.0096 for a/b = 1/200. However, in real 
nanocomposites the transition to the nematic phase is 
hampered by the viscosity of the insulating matrix and 
these systems are inherently out of equilibrium!^ In or- 
der to maintain global isotropy also for > 0i_n, we 
generated oblate spheroid distributions with RSA alone. 
The outcomes are again displayed in Fig. [5] and one can 
appreciate that the difference with the equilibrium re- 
sults for (j) < (/"i-N is quite small and negligible for the 
present aims. 

In analogy to what we have done for the case of pro- 
late objects, it would be useful to find a scaling relation 
permitting to express the </)- dependence of 5c/ D also for 
oblate spheroids, at least for the a/b ^ 1 limit, which is 
the one of practical interest. To this end, it is instruc- 
tive to consider the case of perfectly parallel spheroids 
which can be easily obtained from general result for 
aligned penetrable objects4^ For infinitely thin parallel 
hard disks of radius b one therefore has Ki'L = 2.8/p, 
where FeL = {4/3)'Kb^[l2{5c/ D) + 67r{5c/ Df + 8{5c/ D)^] 
is the excluded volume of the plate plus the penetrable 
shell of critical thickness 5c/2. Assuming that this holds 
true also for hard-core-penetrable-shell oblate spheroids 
with a sufficiently thin hard core, we can then write 



I2{5c/D) + &n{5c/Df + %{5c/Df 



^{b/ay 



(6) 



which implies that 5c/D depends solely on cf>{b/a). As 
shown in the Appendix, where the excluded volume of 
an isotropic orientation of oblate spheroids is reported, 
also the second-order virial approximation gives 5c/ D as 
a function of (f>{a/b) for a/6 ^ 1. Hence, although Eq. © 
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FIG. 5: (Color online) (a): the 5c/D dependence on the volume fraction (j) for impenetrable oblate spheroids with a/b — 1, 1/2, 
1/10, 1/100, and 1/200. Results obtained by RSA alone are also presented, (b): our 5c/D values plotted versus the re-scale 
volume fraction <j)(h/a). The dashed line follows Eq. ([6]) and the solid line follows Eq. ((Tj. Inset: the asymptotic behavior for 
5a/D < 0.1. 



and Eq. (jA13[) are not expected to be quantitatively accu- 
rate, they suggest nevertheless a possible way of rescaling 
the data of Fig. EJa). Indeed, as shown in Fig. E^b), for 
sufficiently high shape anisotropy the data of 5c/ D plot- 
ted as a function of (j){a/b) collapse into a single curve 
(the results for a/b = 1/100 and a/b = 1/200 are com- 
pletely superposed). Compared to Eq. which behaves 
as Sc/D DC [(j){b/a)]~^ for Sc/D <C 1 (dashed line), the re- 
scaled data in the log-log plots of Fig. ^h) still follow a 
straight line in the same range of 6c/ D values but with 
a slightly sharper slope, suggesting a power-law depen- 
dence on (j){a/b). Empirically, Eq. ([5]) does indeed repro- 
duce then the a/b 1 asymptotic behavior by simply 
modifying the small 5c/ D behavior as follows: 



I2a{5c/Df + Q'K{5c/Df + d,{5c/Df 



2.8 



(7) 



where a = 1.54 and (3 = 3/4. When plotted against 
our data, Eq. ([7]) (solid line) provides an accurate ap- 
proximation for 5c/ D in the whole range of (f>{a/b) for 
a/b < 1/100. Moreover, by retaining the dominant con- 
tribution of Eq. ([7|) for 5c/ D < 0.1, we arrive at (inset of 
Fig.EIb): 



C. spheres 

Let us conclude this section by providing an accurate 
expression for 5c/ D also for the case of spherical impen- 
etrable particles. In real homogeneous composites with 
filler shapes assimilable to spheres of diameter in the sub- 
micron range, the cross-over volume fraction (pc is con- 
sistently larger than about 0.1^ so that a formula for 
5c/ D that is useful for real nanosphere composites must 
be accurate in the (j> >0.1 range. For these (j> values the 
scaling relation 5c/ D oc 0"^^^, which stems by assum- 
ing very dilute systems such that 5c/ D ^ 1, is of course 
no longer valid. However, as noticed in Ref. [l5|, the 
ratio 5c/5^f^, where 5nn is the mean minimal distance 
between nearest-neighbours spheres, has a rather weak 
dependence on (j). In particular, we have found that the 
5c data for a/b = 1 in Fig. 2] are well fitted by assuming 
that 5c = 1.65(5nn for > 0.1. An explicit formula can 
then be obtained by using the high density asymptotic 
expression for (5nn as given in Ref. [2^. This leads to: 



5c/D 



1.65(1 



120(2- 



(9) 



which is plotted by the solid line in Fig. Ufa). 



5c/D 



0.15(a/&) 



-,4/3 



(8) 



ANALYTIC DETERMINATION OF THE 
FILLER DEPENDENCIES OF THE 
CONDUCTIVITY 



which applies to all cases of practical interest for plate- 
like filler particles {a/b <C 1 and 5c/D ^ 1). 



With the results of the previous section, we are now 
in a position to provide tunneling conductivity formulas 
of random distributions of prolate, oblate and spherical 
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objects for cr > (Tin, where is the intrinsic conductivity 
of the matrix. Indeed, by substituting Eqs. (|S]), and 
(HI) into Eq. ^ wc obtain 



cr ~ (To exp 



(7 ~ Co cxp < 



(7 ~ (Jo cxp 



2D-f{b/a) 



for prolates, 



(10) 



2D 


'0.15{a/by 




T 


[ <i> \ 


4/3 1 



2D 1.65(1 



e 120(2-0) 



for oblates, (11) 



for spheres. (12) 



From the previously discussed conditions on the validity 
of the asymptotic formulas for ^cjD it follows that the 
above equations will hold when {hja)^ ^ ^ 0.03 for 
prolates, </> > a/& and ajh < 0.1 for oblates, and > 
0.1 for spheres. We note in passing that for the case of 
prolate objects, a relation of more general validity than 
Eq. PH)) can be obtained by substituting Eq. (O into 
Eq. ©. 

Although we are not aware of previous results on a 
for dispersions of oblate (plate-like) particles, there exist 
nevertheless some results for prolate and spherical par- 
ticles in the recent literature. In Ref. [29|, for exam- 
ple, approximate expressions for a for extreme prolate 
{a/b ^ 1) objects and their temperature dependence 
have been obtained by following the critical path method 
employed here. It turns out that the temperature inde- 
pendent contribution to a that was given in Ref. [2^ has 
the same dependence on the particle geometry and den- 
sity of Eq. ([TU)) . but without the numerical coefficients. 
The case of relatively high density spheres has been con- 
sidered in Ref. [3] where In(cr) oc 1/0 has been proposed. 
This implies that 6c/ D oc 1/0, which does not adequately 
fit the numerical results of 6c/ D, while Eq. and con- 
sequently Eq. p2)) . are rather accurate for a wide range 
of values. 

In addition to the 0-dependence of the tunneling con- 
tribution to the conductivity, Eqs. PU | - ([T^ provide also 
estimations for the cross-over value 0c, below which the 
conductivity basically coincides with the conductivity am 
of the insulating matrix. As discussed in Scc lIIIl and as 
illustrated in Fig. [31 0c may be estimated by the value 
such that a ~ cTm, which leads to 



2D -i{b/af 



£, ln{ao/(Jm)' 



for prolate and 

0c ~ 0.15(a/6) 



2D 



^ ln(cro/crm) 



3/4 



(13) 



(14) 



for oblate objects. For the case of spheres, 0c is the root 
of a third-order polynomial equation. Equations (jl3p and 
([T4)) . by construction, display the same dependence on 
the aspect-ratio of the corresponding geometrical perco- 
lation critical densities, as it can be appreciated by com- 
paring them with Eq. Q (prolates) or with the inverse 



of Eq. ([8]) (oblates). However they also show that the 
cross-over point depends on the tunneling decay length 
and on the intrinsic matrix conductivity. This implies 
that if, by some means, one could alter am in a given 
composites without seriously affecting ^ and cto, then a 
change in 0c is to be expected. 



VI. COMPARISON WITH EXPERIMENTAL 
DATA 

In this section we show how the above outlined formal- 
ism may be used to re-interpret the experimental data on 
the conductivity of different nanocomposites that were 
reported in the literature. In Fig. [6l^a) we show the 
measured data of In(CT) versus for polymer compos- 
ites filled with graphene sheets?^ Pd nanosphereS)^ Cu 
nanofiberS)^ and carbon nanotubes4S Equation ^ im- 
plies that the same data can be profitably re-plotted as 
a function of 6c, instead of 0. Indeed, from 



Inia) 



ln((To), 



(15) 



we expect a linear behavior, with a slope —2/^, that is 
independent of the specific value of ao, which allows for a 
direct evaluation of the characteristic tunneling distance 
f. By using the values of D and a/b provided in Refs. [1, 
l47l - [4a | (see also Appendix B) and Eqs. (O, ^ for 
6c, we find indeed an approximate linear dependence on 
6c [Fig- Elb)], from which we extract ^ ~ 9.22 nm for 
graphene, 1.50 nm for the nanospheres, 5.9 nm for the 
nanofibers, and 1.65 nm for the nanotubes. 

We further applied this procedure to several pub- 
lished data on polymer-based composites with nanofibers 
and carbon nanotubes?^ nanospheres;^ and nanoshects 
(graphite and graphene) hence with fillers having 
a/b ranging from ~ 10^'^ up to ^ lO'^. As detailed in 
Appendix B, wc have fitted Eq. ()15|) to the experimen- 
tal data by using our formulas for 6c. The results are 
collected in Fig. [71 showing that most of the so-obtained 
values of the tunneling length ^ are comprised between 

0.1 nm and ~ 10 nm, in accord with the expected 
value rangeJii^'i^i^'^ This is a striking result consid- 
ering the number of factors that make a real composite 
deviate from an idealized model. Most notably, fillers 
may have non-uniform size, aspect-ratio, and geometry, 
they may be oriented, bent and/or coiled, and interac- 
tions with the polymer may lead to agglomeration, segre- 
gation, and sedimentation. Furthermore, composite pro- 
cessing can alter the properties of the pristine fillers, e.g. 
nanotube or nanofibcr breaking (which may explain the 
downward drift of ^ for high aspect-ratios in Fig. [7]) or 
graphite nanosheet exfoliation (which may explain the 
upward shift of ^ for the graphite data). In principle, 
deviations from ideality can be included in the present 
formalism by evaluating their effect on 6c^ It is how- 
ever interesting to notice that all these factors have of- 
ten competing effects in raising or lowering the composite 
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FIG. 6: (Color online) (a): natural logarithm of the conductivity a as a function of the volume fraction (f) for different polymer 
nanocomposites: Graphene-polystyrene,— Pd nanospheres-polystyrene,— Cu nanofibers-polystyrene,— and single-wall carbon 
nanotubes-epoxy.— When, for a given concentration, more then one value of a was given (as in Refs. [Uliil), the average of 
In(cr) was considered, (b), the same data of (a) re-plotted as function of the corresponding critical distance 5c- Solid lines are 
fits to Eq. (fT5|) . 
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FIG. 7: (Color online) Characteristic tunneling distance ^ values for different polymer nanocomposites as extracted form 
Eq. (|15p applied to the data of Refs. [47II5H (low structured carbon black and metallic nanosphere composites), Refs. [i^-lsol] 
(nanofiber and carbon nanotube composites), and Refs. P lsolls^ (nanographite and graphene composites). 



conductivity, and Fig.lHlsuggcsts that on the average they 
compensate each other to some extent, allowing tunnel- 
ing conduction to strongly emerge from the (^-dependence 
of (T as a visible characteristic of nanocomposites. 



VII. DISCUSSION AND CONCLUSIONS 

As discussed in the introduction, the theory of con- 
ductivity in nanocomposites presented in the previous 
sections is based on the observation that a microscopic 



mechanism of interparticle conduction based on tunnel- 
ing is not characterized by any sharp cut-off, so that 
the composite conductivity is not expected to follow 
the percolation-like behavior of Eq. ([T]). Nevertheless, 
we have demonstrated that concepts and quantities per- 
tinent to percolation theory, like the critical path ap- 
proximation and the associated critical path distance 
(5c, are very effective in describing tunneling conduc- 
tivity in composite materials. In particular, we have 
shown that the (geometrical) connectivity problem of 
semi-penetrable objects in the continuum, as discussed 
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FIG. 8: (Color online) Conductance versus dependence for 
a carbon black-quaternized Poly(4-vinylpyridine) composite 
for different humidities. Adapted from Ref. [57| 



in Sec. IIVI is of fundamental importance for the under- 
standing of the filler dependencies (0, D, and a/6) of cr, 
and that it gives the possibility to formulate analytically 
such dependencies, at least for some asymptotic regimes. 
In this respect, the body of work which can be found on 
the connectivity problem in the literature finds a straight- 
forward applicability in the present context of transport 
in nanocomposites. For example, it is not uncommon to 
find studies on the connectivity of semi-penetrable ob- 
jects in the continuum where the thickness (5/2 of the 
penetrable shell is phenomenologically interpreted as a 
distance of the order of the tunneling length ^, 38i39,55^56 
This interpretation is replaced here by Eq. ([3]) which pro- 
vides a clear recipe for the correct use, in the context of 
tunneling, of the connectivity problem through the crit- 
ical thickness (5c/2. Furthermore, Eq. ([3]) could be ap- 
plied to nanocomposite systems where, in addition to the 
hard-core repulsion between the impenetrable particles, 
effective inter-particle interactions are important, such as 
those arising from depletion interaction in polymer-base 
composites. In this respect, recent theoretical results on 
the connectivity of polymer-nanotube composites may 
find a broader applicability in the present context 

It is also worth noticing that, although our results 
on the filler dependencies of 6^ for prolate objects with 
a/& ^ 1 can be understood from the consideration of 
excluded volume effects (e.g. second virial approxima- 
tion), the corresponding be formulas for the oblate and 
spherical cases are empirical, albeit rather accurate with 
respect to our Monte Carlo results. It would be there- 
fore interesting to find microscopic justifications to our 
results, especially for the case of oblates with ajh <Si 1, 
which appear to display a power-law dependence of Sc on 
the volume fraction [Eq. 

Let us now turn to discuss some consequences of the 



theory presented here. As shown in Sec.|V]the cross-over 
volume fraction (f>c depends explicitly on the conductivity 
CTm of the insulating medium, leading to the possibility 
of shifting (j)c by altering (Jm. Formulas P3)) and 
were obtained by assuming that the transport mecha- 
nism leading to was independent of the concentration 
of the conducting fillers, as it is the case for polymer- 
based nanocomposites, where the conduction within the 
polymer is due to ion mobility. In that change in 

CTm, and so a change in 0c, could be induced by a change 
in the ion concentration. This is nicely illustrated by 
an example where a conductive polymer composite with 
large ionic conductivity was studied as a material for hu- 
midity sensors 1^ This consisted of carbon black dispersed 
in a Poly(4-vinylpyridine) matrix which was quaternized 
in order to obtain a polyclcctrolyte. Since the absorbed 
water molecules interact with the polyelectrolyte and fa- 
cilitate the ionic dissociation, higher humidity implies a 
larger ionic conductivity. In Fig. [5] we have redrawn Fig. 
4 of Ref. [131 terms of the conductance as a function of 
carbon black content for different humidity levels. Con- 
sistently with our assumptions, one can see that with 
the increase of humidity the matrix intrinsic conductiv- 
ity is indeed shifted upward, while this has a weaker effect 
on the conductivity for higher contents of carbon black, 
where transport is governed by interparticle tunneling (a 
slight downshift in this region is attributed to enhanced 
interparticle distances due to water absorption) . The net 
effect illustrated in Fig. [5] is thus a shift of the crossover 
point (pc towards higher values of carbon black content. 
It is worth noticing that the explanation proposed by the 
authors of Ref. [57[ in order to account for their finding 
is equivalent to the global tunneling network/crossover 
scenario. 

Another feature which should be expected by the 
global tunneling network model concerns the response of 
the conductivity to an applied strain s. Indeed, by using 
Eq. ([3]), the piezoresistive response F, that is the relative 
change of the resistivity a^^ upon an applied e, reduces 
to: 



^ ^ c;in(g ^) _ j^/0£\ dlnjSc 
de \ a / de 



(16) 



In the above expression d\n{Sc)/de = 1 for fillers hav- 
ing the same elastic properties of the insulating matrix. 
In contrast, for elastically rigid fillers this term can be 
rewritten as [d\n{5c) / d\B.{(t))]dln.{(j)) / de , which is also ap- 
proximatively a constant due to the 5c dependence on 
as given m Eqs. (gl), (HI), and and to d\n{(j))/de ~ -1. 
Hence, the expected dominant dependence of F is of the 
form F oc ln(l/(T), which has been observed indeed in 
Refs.Hilll. 

Finally, before concluding, we would like to point out 
that, with the theory presented in this paper, both the 
low temperature and the filler dependencies of nanocom- 
posites in the dielectric regime have a unified theoretical 
framework. Indeed, by taking into account particle ex- 
citation energies, Eq. ^ can be generalized to include 
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inter-particle electronic interactions, leading, within the 
critical path approximation, to a critical distance 6c 
which depends also on such interactions and on the tem- 
perature. The resulting generalized theory would be 
equivalent then to the hopping transport theory cor- 
rected by the excluded volume effects of the impenetra- 
ble cores of the conducting particles. An example of this 
generalization for the case of nanotubc composites is the 
work of Ref. 

In summary, we have considered the tunneling- 
percolation problem in the so far unstudied intermediate 
regime between the percolation-like and the hopping-like 
regimes by extending the critical path analysis to systems 
and properties that are pertinent to nanocomposites. 
We have analyzed published conductivity data for sev- 
eral nanotubes, nanofibers, nanosheets, and nanospheres 
composites and extracted the corresponding values of the 
tunneling decay length ^. Remarkably, most of the ex- 
tracted ^ values fall within its expected range, showing 
that tunneling is a manifested characteristic of the con- 
ductivity of nanocomposites. Our formalism can be used 
to tailor the electrical properties of real composites, and 
can be generalized to include different filler shapes, filler 
size and/or aspect-ratio polydispersity, and interactions 
with the insulating matrix. 
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Appendix A: Excluded volumes of spheroids and 
spherocylinders with isotropic orientation 
distribution 

The work of Isiharai^ enables to derive closed relations 
for the excluded volume of two spheroids with a shell of 
constant thickness and for an isotropic distribution of the 
mutual orientation of the spheroid symmetry axes. Given 
two spheroids with polar semi-axis a and equatorial semi- 
axis b, their eccentricity e are defined as follows: 



1 - 



&2 



e=\ l- 



62 



for prolates, (Al) 
for oblates. (A2) 



If the mutual orientation of the spheroid symmetry axes 
is isotropic, the averaged excluded volume of the two 



spheroids is then (valid also for more general identical 
ovaloids) 



2V 



MF 



(A3) 



where V is the spheroid volume and M and F are two 
quantities defined asi^ 



AI =2-Ka 



l + ii^lnfi±^ 

2e \l-e 

arcsin e 



F =27ra6^Vl - e2 



for the case of prolate {a/b > 1) spheroids and 



M ^2nb( Vl - e 



F =2tt¥ 



1 



(1 



2e 



In 



1 - e 



(A4) 
(A5) 

(A6) 
(A7) 



for the case of oblate (a/b < 1) spheroids. 

If now the spheroids are coated with a shell of uniform 
thickness d (d = S/2), then the averaged excluded volume 
of the spheroids plus shell has again the form of (jA3|: 



tot 

exc 



2Vd 



2ti ' 



(A8) 



and by constructing the quantities Vd^ Md , and Fd from 
their definition in Ref. [601 (^^^ ^ei- [2l| for a similar 
calculation), one obtains: 



tot 
cxc 



AdF 



+ 8dHl + —d^. (A9) 

TT 3 



In the cases of extreme prolate {a/b ^ 1 and 6/a ^ 1) 
and oblate {a/b ^ 1 and S/b <^ 1) spheroids, the total 
excluded volume reduces therefore to 

Kxc = K=xc + 2na^S, (AlO) 

for prolates and 

Kxc = K=xc + 47r52,5 + Tr^b^S/2 + 2ti%6^ , (All) 

for oblates. Within the second-order virial approxi- 
mation, the critical distance 5c is related to the vol- 
ume fraction <f> through <j) ~ l//AT4xc, where AT4xc = 
Kxc ~ Ksxc- From the above expressions one has then 
[D = 2max(a,fe)]: 



for prolates and 



35c/ D' 
(4/3)(a/6) 



{8 + tt'^)5c/D + 8t:{5c/DY' 



(A12) 



(A13) 



for oblates. 
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For comparison, we provide below the excluded vol- 
umes of randomly oriented spherocylinders. These are 
formed by cylinders of radius R and length L, capped 
by hemispheres of radius R. Their volume is ^ = 
(4/3)7ri?3[i + (3/4)(L/i?)]. The excluded volume for 
spherocylinders with isotropic orientation distribution 
was calculated in Ref. 1401 and reads 



Vexc = 



l + l{L/R) + ^iL/Rr 



(A14) 



For the high aspect-ratio limit {L/R^ 1), when d/L <C 
1, the total excluded volume minus the excluded volume 
of the impenetrable core is 



(A16) 



which coincides with the last term of Eq. (|A10p if d = (5/2 
and a = i? + L/2~ L/2. Furthermore, the second-order 
virial approximation ~ V/Ayoxc) gives: 



The excluded volume V^'^^' of spherocylinders with a shell 
of constant thickness d = 5/2 is hen: 



V^^l = ^^{R + df 



1 + ^ 



L 



4:\R + d 



L 



32 \R + d 



(A15) 



(A17) 



which has a numerical coefficient different from Eq. (|A12[) 
because for spherocylinders V ~ 2'iTab^ (for a/h ^ 1) 
while for spheroids V = 47ra6^/3. 
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Appendix B: Supplementary material. Conductivity versus critical distance plots 

We show in the following Figs. IQIITTI the complete set of plots of the natural logarithm of the sample conductivity 
(T as a function of the geometrical percolation critical distance for different polymer nanocomposites, as used to 
obtain the ^ values of Fig. 7 of the main article. In collecting the published results of a versus (j), we have considered 
only those works where a/b and D = 2max(a, 6) were explicitly reported. In the cases of documented variations of 
these quantities, we used their arithmetic mean. The 4> dependence of the original published data was then converted 
into a Sc dependence as follows. For fibrous systems (nanofibers, nanotubes), the filler shape was assimilated to 
spherocylinders, while for nanosheet systems it was assimilated to oblate spheroids. For prolate fillers 6c was obtained 
from Eq. (5), for oblate fillers the values for 5c were obtained from Eq. (8), while for spherical fillers, the values of 5c 
were obtained from Eq. (9). 

Since the model introduced in the main text is expected to be representative only if (p is sufficiently above <f>c to 
consider the effect of the insulating matrix negligible, for a given experimental curve, higher (j) data were privileged, 
and lower density points sometimes omitted when deviating consistently from the main trend. The converted data 
were fitted to Eq. (15) of the main text and the results of the fit are reported in Figs. CTTTlbv sohd hues. The results 
for 2/^ and In(ao) are also reported in the figures. As it may be appreciated from Figs. [MT71 in many instances the 
experimental data follow nicely a straight line, as predicted by Eq. (15) of the main text, while in others the data are 
rather scattered or deviate from linearity. In these latter cases, the fit to Eq. (15) is meant to capture the main linear 
trend of In(tT) as a function of 5c. It should also be noticed that, in spite of the rather narrow distribution of the 
extracted ^ values reported in Fig. 7 of the main text, the values of the prefactor ctq obtained from the fits are widely 
dispersed. This is of course due to the fact that, besides intrinsic variations of the tunneling prefactor conductance 
for different composites, interpolating the data to 5c = leads to a large variance of co even for minute changes of 
the slope. We did not notice any significant correlation between the extracted ^ and ctq values. 
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FIG. 12: plots of Intj as a function of 5 for different polymer-nanotube composites (cont.) 
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FIG. 13: plots of In cr as a function of 5 for different polymer-nanotube composites (cont.) 
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FIG. 14: plots of In cr as a function of 5 for different polymer-nanotube composites (cont.) 



21 



Flandin et al., J. Appl. Polym. Sci. 76, 894 (2000) 
CB D=300 nm 




S [nm] 

Rubin et al., Phys. Rev. B 59, 12196 (1997) 
CB D=320 nm 








-14 I , , , , , 1 

2 4 6 8 10 12 

5 [nm] 

Mohanraj et al., Polym. Compos. (2007) DOI 10.1002/pc.20236 
Ni-Cu D=70 nm 



Nakamura et al., Jpn. J. Appl. Phys., Part I 36, 5163 (1997) 




5 10 15 



5 [nm] 



Kubat et al., Synthetic Met. 54, 187 (1993) 
Pd D=300 nm 

-5 I 




-25 I ' ' ' ' ' ' 1 

15 17 19 21 23 25 27 29 

5 [nm] 

Untereker et al., Appl. Mater. Interfaces 1, 97 (2009) 
Pt D=50 nm 





5 


E 





o 


-5 


CO 






-10 




-15 




-20 




-25 




-30 



50 



y = -0.3808X + 28.001 

70 90 110 

5 [nm] 



130 



150 




FIG. 15: plots of In £7 as a function of S for different polymer-nanospheres composites 
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FIG. 17: plots of Incr as a function of 5 for different polymer-nanosheet composites (cont.) 



